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Abstract 

We consider the problem of estimating the 2D vector displacement field in a 
heterogeneous elastic solid deforming under plane stress conditions. The problem 
is motivated by applications in quasistatic elastography. From precise and accu¬ 
rate measurements of one component of the 2D vector displacement field and very 
limited information of the second component, the method reconstructs the second 
component quite accurately. No a priori knowledge of the heterogeneous distribu¬ 
tion of material properties is required. This method relies on using a special form 
of the momentum equations to filter ultrasound displacement measurements to pro¬ 
duce more precise estimates. We verify the method with applications to simulated 
displacement data. We validate the method with applications to displacement data 
measured from a tissue mimicking phantom, and in-vivo data; significant improve¬ 
ments are noticed in the filtered displacements recovered from all the tests. In 
verification studies, error in lateral displacement estimates decreased from about 
50% to about 2%, and strain error decreased from more than 250% to below 2%. 


1 Introduction 


Ultrasound elasticity imaging (UEI), or elastography, is a rapidly growing field, pri- 
marily due to the increasing interest in the non-invasive quantification of the mechan¬ 
ical properties of soft tissues |Gao et ah, 1996[ |Ophir et ah, 19991 |Parker et ah, 2005 


Greenleaf et ah, 2003[ |Parker et ah, 20lT| |Doyley, 2012 J. A crucial step in UEI is the 
estimation of a tissue’s motion. This motion is typically generated by a quasistatic 
compression, and the deformation of the tissue is measured with ultrasound. The mea¬ 
sured deformation and an appropriate mathematical model can then be used to infer 
the mechanical properties of the soft tissue. 

Ultrasound measurements of tissue deformation are typically much more precise in 
the “axial” direction (i.e. the direction of the ultrasound beam), than in the so-called 
“lateral” direction, the image direction transverse to the direction of the ultrasound 
beam. The origin of the difference in precision is due to the anisotropy of the ultrasound 
point spread function (PSF). The width of the peak of the radio frequency (RF) PSF, in 
the axial direction (the direction of sound propagation) is roughly an order of magnitude 
smaller than that in the lateral direction |Szabo, 2004). 
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With only one component of the displacement field measured precisely, it is im¬ 
possible to compute precise estimates of the full strain tensor field, which would oth¬ 
erwise be useful information to researchers who do strain imaging. Furthermore, one 
goal of measuring the displacement field is as input to an inverse problem to estimate 
material property distributions |Ba rbone and Qberai, 2010] , |Doyley, 2012| . Computing 
the material properties from only one component of the displacement field is compu¬ 
tationally intensive |Ba rbone and Qberai, 2010] , |Doyley, 2 012 . In that situation, an 
iterative inversion approach is required, which may takes hundreds to thousands of 
iterations to converge. Such an iterative inversion code usually needs displacement 
boundary conditions to drive the forward problem, and using the noisy measurements 


as boundary conditions can corrupt the material property estimates Dord et ah, 2012 


Goenezen et ah, 2012[ Richards et ah, 2 009 . On the other hand, efficient direct meth¬ 


ods to compute material properties from full vector displacement field measurements 
exist |Albocher et ah, 2009[ |Barbone et ah, 2010[ |Albocher et ah, 2014| . The present 


contribution uses an inverse problem approach that aims to improve lateral displace¬ 
ment estimates from ultrasound measurements. 

Much research has been performed in order to recover better lateral displacements 
both in the ultrasound elasticity imaging and the ultrasound blood flow imaging field. 


Geiman et ah, 2000, Konofagou and Ophir, 2000, 

Chen et ah, 2004| Ebbini, 2006, Brusseau et ah, 2008 

Deprez et ah, 2009, Kim et ah, 2011, Rivaz et ah, 2011, 

Mailloux et ah, 1987, Basarab et ah, 2008 


the authors use various forms of interpolation within standard or modified block-matching 
algorithms to estimate subsample lateral displacements. These approaches try to make 
the best use of available data as is, but do not overcome the fundamental difficulty posed 
by the anisotropic PSF. An alternative approach is to use beam steering to measure dis¬ 
placements at arbitrary angles, and then use the arbitrary angled displacement estimates 
to determine the axial and lateral displacements as described in [Tanter et ah, 2002 


|Techavipoo et ah, 2004[ |Rao et ah, 2007[ |Hansen et ah, 20T0] , |Xu and Varghese, 2013[ 


Abe ysekera et ah, 2012 [. This method uses additional image data to overcome the lim¬ 


itations of a single PSF. A third approach, pursued by [Korukonda and Doy ley, 2011 


Korukonda and Doyley, 2012[ |Jensen and Munk, 1998[ |Liebgott et ah, 2007| , is to use 


novel beam forming techniques to modulate the RF signal in the lateral direction so 
that it contains more phase information. This approach tries to reduce the anisotropy 
of the PSF. A fourth approach, pursued by [Ri chards et ah, 2009[ [Lubins ki et ah, 1996] , 
|Skovoroda et ah, 1998[[O’Donnell et ah, 2001| , is to constrain the displacements to sat¬ 
isfy the incompressibility condition. Many different types of soft tissue can be assumed 
to be incompressible because they are mostly composed of water. 

Among the existing approaches, the incompressibility processing method introduced 
in [Lu binski et ah, 1996| is arguably the most similar in spirit to the method presented 
here to reconstruct better lateral displacements. In that paper, the authors assumed 
that the deformation is linear, incompressible, and plane strain in character. With these 
assumptions, they formulated a minimization problem where the objective was to find 
lateral displacements that minimized the error between the measurements and predicted 
lateral displacements that satisfied the modeling assumptions. The incompressibility 
processing described in [O’Donn ell et a h, 2001 is basically the same as the method de¬ 
veloped in [Lubinski ej^al., 1996], with the addition of weighting functions to the lateral 
displacement reconstruction equations. These weighting functions are based on the abso¬ 
lute value and the gradient of the cross-correlation function. In [Skovoroda et ah, 1998j , 
the method developed in [Lubinski et ah ,JL99 61 is extended to accommodate large plane 
strain deformations. Finally, in [Richards et ah, 2009| , the incompressiblility constraint 
is used to improve the lateral and elevational displacements for a 3D deformation. 
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In this paper, we introduce a fundamentally new approach to evaluating lateral 
displacements. We will be using a spatial regularization term within an inverse problem 
formulation to adaptively enforce, and locally relax, a special form of the momentum 
conservation equation on the measured displacement field. Special weighting functions 
will be used to place more emphasis on the axial component of the displacement field. 

The rest of the paper is arranged as follows: section [2] describes the new formulation 
to estimate better lateral displacements; we refer to this method as “SPREME,” which 
stands for SParse RElaxation of the Momentum Equation. Section [3] describes the 
simulation experiments used to test and verify the displacement filtering algorithm. 
Section [4] describes the phantom experiments used to validate the performance of the 
displacement filtering algorithm on ultrasound measured data. Section [5] describes the 
in-vivo experiments used to further test the performance of the displacement filtering 
algorithm on patient data. Finally, a brief summary and conclusion is given in section 
[6] along with possible future directions of the work. The appendix describes a simpler 
formulation that produces apparently improved displacement and strain fields. Despite 
appearances, however, the resulting strain and displacement fields are non-physical. 


2 Formulation 


We suppose we are given measured displacements in a 2D region D. We further suppose 
the measurements of the axial displacement, u y , are significantly more precise and accu¬ 
rate than the measurements of the lateral displacement. We make use of the assumption 
that within the observed plane, the body’s deformation may be well approximated by the 
plane stress approximation. These assumptions lead to the following problem statement: 


Problem 2.1. Given 2D measured displacement u m {x,y), for all x G D, find the 2D 
displacement vector u{pc,y), and the 2D linearized strain tensor e(x,y) such that: 



Ko[u) = l||(Wm ~u)\\ 2 t 


— V / (y>m R') ’ T(l6 m u ) dD 
z Jn 

is minimized, and: 

e = V s u 

and: 

V • A = 0 a.e. in D 

where: 

A(e) = 2tr(e)I + 2e. 

and: 

V s u = i (v« + (Vtt) T ) 


With diagonal T, equation 0 may be written explicitly as: 


O \^ xx ^ U (rn)x u x) + r ^yy{ u (m)y u y) ] d£l. 
J O 


(1) 

( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 


In equation ([6]), the weights (T xx , T yy ) allow more importance to be placed on the ac¬ 
curate component of the displacement field while calculating the predicted displacement 
field. In equation ([3]), a.e. means almost everywhere , and implies that the constraint is 
enforced everywhere in the domain except on a set of measure zero (i.e. along curves 
and at isolated points). 
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To motivate constraint ft we consider a thin sheet in the (x, y) plane made of solid 
material that is linearly elastic, isotropic, incompressible, and simply connected. We 
further assume that the deformation of the material is small, quasistatic, and approx¬ 
imately plane stress in character. For these modeling assumptions, linear strain e and 


Cauchy stress a are defined and related as: 

e ij = \(di u j + dj u i) (*,.7 = 1,2) (7) 

er = —pi + 2 pe (8) 

where 0 t = For a plane stress deformation: 

0zz = ~P + 2ye zz = 0. (9) 

For an incompressible material: 

^xx T £yy T ^zz ~ 0- (10) 

Substituting (10) into 0 gives: 

P = Qfl^Exx T €yy)- (H) 

Using in ([§]) yields: 

<j — ‘2[i{e xx T T 2/ie (12) 

fiA. (13) 


The equilibrium equation, with no body force, is V • cr 


0, or using (13): 


V • (/iA) = 0 in Q. 


(14) 


If p(x,y) is piecewise constant, then V • A = 0 a.e. in Q. 


2.1 Uniqueness of solution 

Given an axial displacement measurement u y , the conditions necessary to guarantee 
the uniqueness of the lateral displacements computed using the equilibrium equation 
constraints will be derived in this section. This derivation begins with the assumption 
that the material has a piecewise homogeneous shear modulus distribution. This means 
that the material’s domain, fl, can be divided into n different subdomains, and the 
shear modulus in each subdomain is constant. The shear modulus is not constant in the 
interface between the subdomains, however. An illustration of this type of domain, for a 
special case when n — 2, is shown in figure 0- The conditions necessary for uniqueness 
of the lateral displacements in each subdomain will be derived next. 


2.1.1 Solution in a single subdomain 

Lemma 2.1. Given u y (x,y)in subdomain fl( n \ in which fi = constant. Then the 
general solution of equation Q) - (^) for u x {x,y), with (x,y) E is: 

u x (x,y ) = u p x (x,y) + c^x 2 -4y 2 ) + c 2 x + c 3 y + c 4 (15) 

Proof. Since // = constant in Q (n l, we may expand equation (pj) in terms of displacements 
via equations 0 and 0, which results in the following equations: 

^ u x,xx ^x,yy — ^^y,yx 

3 ^x,xy — ^y^xx T 4 u y,yy 
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(16) 

(17) 




All the terms on the right hand side of both these equations are supposed to be 
known, since by assumption here, u y is known precisely. The terms on the left hand 
side involving the lateral displacements, u x , are unknowns. We now show that these 
equations lead to a solution for u x that has at most four unknown constants. We begin 
by splitting the total solution into homogeneous and particular parts. To that end, we 
let: 

(18) 


u x = vL + u‘ 


Here u x is any particular solution of (16) and s which, by assumption, has no free 
constants or undetermined functions. If they appear in the process of solving (16) and 


then they shall be set to arbitrary values (e.g. zero) without loss of generality. Any 
part of the general solution (18) that is undetermined by u y is thus contained in the 


homogeneous solution of (16) and ©• The homogeneous solution satisfies: 

= 0 


Ay h -U 7A 

1 x,yy 


3 U\ 


x,xy 


= 0 . 


Integrating (20) twice yields: 


Substituting (21) into (19) gives: 


u x = f(x)+g(y). 


4 f"(x) + g"(y) = 0. 


Equation (22) implies: 


4/" 0*0 = —g"{y) = 8cj. 


Integrating (23) for f{x), and g(y) yields: 

f{x) — c 1 x 2 + c 2 x + c' 4 
g(y ) = -4c, y 2 +c 3 y + c' 


Combining (24), (25), in (21) gives: 

U X = C 1 (% 2 - 4 y 2 ) + C 2 x + c 3 y + c 4 , 

where c, = c' + c '. 


Therefore, the total solution is of the form (15). 


(19) 

( 20 ) 

( 21 ) 

( 22 ) 

(23) 

(24) 

(25) 

(26) 

□ 


Equation (15) shows that within any homogeneous subregion of the body, given 


the axial displacement u y (x,y ), then the lateral displacement u x (x,y ) has at most 4 
undetermined constants. 


2.1.2 Solution in connected subdomains covering entire domain of interest 

Interestingly, the form of the homogeneous solution is completely independent of the 
given u y (x,y). Thus we may conclude almost immediately: 

Lemma 2.2. Given u y (x,y ) in Q = in which /i(x,y) is piecewise constant, such 

that 

M(z,y) = {/i (n) if(x,y)en^. (27) 

Then the general global solution of equations |1|) - |5|) for u x (x,y) with (x,y) € = 

UOM is: 

u x (x, y) = u p x (x, y) + c x ( x 2 - 4y 2 ) + c 2 x + c 3 y + c 4 (28) 

Here, u x (x,y) is determined entirely by u y {x,y), while c 1 ,...c 4 are undetermined con¬ 
stants. 
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Proof. We consider several contiguous subregions, each of which is homogeneous. From 


lemma 2.1, the solution within subregion n is given by: 

u^Xx, y ) = *(a, y ) + c[ n >(* 2 - Ay 2 ) + c^x + c^y + c< n >. (29) 

If we consider the displacement u x (x,y) to be continuous at a non-special interfac^] 
between subregion n and adjacent subregion m, and we assume that u^ P = u^ p is 
continuous on the interface, then we conclude: 


c (n) _ c (m)_ 


(30) 


To see this, note that and u^ P — u^ p leads to the equation: 


[x 2 — Ay 2 x y l] 




(31) 


with c* = c — c( m ). We note that (31) holds at every point on the interface. Provided 


that the functions x 1 — Ay , x , y , 1 are linearly independent on the interfaces, then the 


only solution of (31) is c* = 0. Interfaces on which these functions are linearly dependent 
are special cases, and we call them “special interfaces.” 

Since the Cj’s take the same value in every subregion, we may consider the homoge¬ 
neous functions to be defined globally. Thus, the global solution is: 


u x {x,y) = u p x (x,y) + c^x 2 -Ay 2 ) + c 2 x + c 3 y + c 4 


(32) 


□ 


Lemma 2.2 and equation (32) implies that by knowing the axial displacement through¬ 


out the domain in a piecewise homogeneous material, the lateral displacement may be 
determined up to four coefficients. Hence, if sufficient information exists in the mea¬ 
surements to determine these four coefficients accurately, then the lateral displacement 
field may be determined accurately throughout the domain. 

In the next section, we derive a variational formulation that permits this estimation 
in a natural way. 


2.2 SPREME variational formulation 

In practice, we shall enforce constraint ^ approximately as a penalty, in the form of a 
regularization term. To determine the appropriate form of that term, we refer to Figure 
0 and consider a continuous y of the form: 

{ yi x E Oi 

fi 2 x G (33) 

continuous monotone x G 

where Pt c is the set of all points within | from curve C. We assume there is a constant 
C\ so that |V/i| < Xt f° r a U x £ 

2 A “special interface” would be one for which there exists [01,02,03,04] = [0*^2,03,04] ^ [0,0,0, 0] 
s.t. c*(x 2 — Ay 2 ) + c\x + c*y + c* = 0 for all (x, y ) on the interface. 
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Figure 1: A piecewise continuous distribution of elastic modulus represented as the 
limit of a continuous distribution. The function y(x,y) has the constant value /ii for 
(x,y) G fix, has the constant value y 2 for ( x,y ) G and is continuous and monotone 
in a narrow intermediate region, Q c . 

Suppose V • (yA) = 0 everywhere in fi, then: 


{ /iiV • A = 0 x G fii 

y 2 V • A = 0 x G fi 2 • 

/xV • A A V y = 0 x G fic 

In order to form the penalty term, we assume (l34j) is satisfied and consider: 


(34) 


[ |V • A| Q dQ = [ 

Jn Jq. ( 


V/x 
’ M 


dQ r < 


/ 


iAnv^i c 


■dQ c < 


I 

jQr 


\ cx 

E 2 


“w IT 




Y ) dn c (35) 
(36) 


where // mm = min{/ii,/i 2 } is the minimum value of /i in fi c , Ch = (C'i | A|)/ /J mm , |A| is 
the largest eigenvalue of A, and is the perimeter of curve C; see Figure (JTJ) . 

We see in (36) that for piecewise constant y (i.e. when h —>► 0), 


/1 


IV-Al“dfi = 0 


for 


0 < a < 1. 


(37) 


Equation (37) shows that we may enforce the constraint V • A = 0 a.e. by choosing to 


enforce the integral condition (37) with 0 < a < 1. 


The derivative of the absolute value function is undefined at the origin, so in practice 
we introduce a small positive constant, 5, and approximate (|37|) as: 


[ |V-A(e)| a dtt^ [ 

Jq Jq 


(V-A (e k )Y 


(V-A(e (fc_1) )) 2 + <5 


dQ 


where: 


1 a 

n— 1- 

2 


(38) 


(39) 


Note that for 0 < a < 1 


1 > n > k. 


In Equation (38), k denotes an iteration counter. As k approaches infinity, we expect 

(k— 1) j 

e to approach e . In this limit, for 6 = 0, the approximation in (38) is then exact. 


The S in the denominator is a small numerical constant that prevents the denominator 
from being zero when V • A = 0. 

Using (38) instead of <§> we re-formulate the problem as: 
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Problem 2.2. Given u m (x,y), x G Gt and e^ k l \x,y), find u k , and e k that minimizes: 


where: 


and: 


and: 


and: 


7R 


7l[U, e\ = 7T 0 + 7T C + 7 T r 


= \ [ ( U m ~ U k ) ■ T (U m - U k ) (IQ 

z J Q 


^ [ (e k -V s u k ) 2 dn 
Jq 


= P 

' c 2 


7Tr 2 


lL a{x) ( v ' A< 


• A(e k )) dn 


a{x) = 


OL n 


(V-A(e (fc_1) )) +5 


(40) 

(41) 

(42) 

(43) 

(44) 


This yields a quadratic minimization problem at each iteration. We initialize itera¬ 
tions using e® = 0. Equation (40) contains a displacement data matching term ( 7 r G ), 
a compatibility term ( 7 r c ), and a regularization term ( 71 ^). 7r 0 minimizes the mismatch 
between the measured displacement field and the filtered displacement field, n c con¬ 
strains the calculated strains to satisfy the compatibility equation. (3 is a constant that 
is used to regulate how strongly the compatibility condition is enforced. n R is used 
to constrain the strains to satisfy the equilibrium equations almost everywhere in the 
domain. a(x) maybe interpreted as a spatially varying regularization constant that 
regulates how strongly the equilibrium equation is enforced locally. 

The way the regularization, 7 r H , works can be best explained with figure ([lj. In 
this figure, where fi is constant (i.e. in the background (/i 2 ) and in the inclusion (/xi)), 
V • A(e ) will be small meaning that a(x) will be large, so the equilibrium equation 
will be strongly enforced in these regions. Therefore, u x can be confidently determined 
there up to the homogeneous solution. Where ji is not constant (i.e. in a thin 
region along the interface between the inclusion and the background), V • A(e ) 
will be large, meaning that a(x) will be small, so the equilibrium equation will not be 
strongly enforced. Thus the strains are permitted to change rapidly there. Continuity 
of the lateral displacements is enforced there through the n c term. This “spatially 
adaptive” regularization thus allows us to enforce the equilibrium equation everywhere 
in the domain, except near material interfaces. 

We note that if the goal is to evaluate the strain tensor only, then it is tempting 
to reformulate the problem without the displacement variable, u. Such a formulation 
is presented in appendix & There, it is shown that the strains obtained from the 
formulation are not compatible, and hence not physical, for some simulation experiments. 

Equation (40) is the SPREME functional. Our implementation minimizes 7 r[te, e] by 
the finite element method (FEM) discretization. The weak form is derived first by 
computing the following directional derivatives: 


D u n:v 


D e n : w 


d 

dc 


7 t[u + cv, e] = 0 

c=0 


d 

dc 


7 t[u, e + cw] — 0 


(45) 

(46) 


c=0 









where v and w are admissible variations of u and e, respectively, and e = e k . Taking 
the directional derivatives yields the following system of equations: 



v dVt — (3 


In 


V s u] (Vv) dfl = 0 


for all v 


(47) 


f3 f [e — \7 s u} w dft + f <a(a?){V • A(m)} • {V • A(e)} = 0 for all w (48) 

Jn Jn 


Equations (47) and (48) were discretized with bilinear shape functions to yield a 
system of linear equations. A new finite element within an in-house FEM code was 
created to solve the system of linear equations following standard FE methods e.g. 
|Hughes, 2000| . 


3 Verification of computational implementation 

The purpose of verification is to confirm that, given data that satisfies the modeling 
assumptions, a computational implementation produces the correct solution that is con¬ 
sistent with the assumed model. In this section, therefore, the results of two different 
simulation experiments conducted to test the implementation of the SPREME formulation 
will be described. For the two experiments performed, reference analytical displacement 
and strain fields are available, and they were used to evaluate the quality of reconstructed 
displacements and strains from the SPREME formulation. 

3.1 Reference analytical solution 

An exact analytical solution of the elasticity equations is used to generate reference 
displacement data. The solution corresponds to that of an infinite elastic sheet with a 
perfectly bonded circular inclusion. The sheet is subject to uniform stress at infinity. 
The problem was nondimensionalized using the region of interest size as the reference 
length, and the background shear modulus as the reference stress. The shear modulus of 
the inclusion and background were (i T — 3, and fi B = 1, respectively. The Poisson’s ratio 
was taken to be v — 0.5, consistent with the incompressibility assumption. The radius 
of the inclusion was a — 0.25. The solution to the elasticity equation was obtained using 
the methods described in |Honein, 1990| . 

3.2 Uniaxial tension test 

For the first experiment the sheet was subject to uniaxial tension in the y direction as 
depicted in figure ([ 2 ]). The stress at infinity was taken to be a yy = <j°° = 6 x 10 -2 . 
The displacement field corresponding to this deformation was evaluated in a subregion 
around the inclusion. This subregion was discretized into 50 x 50 square elements in 
the x and y directions, and the axial and lateral displacements were evaluated at each 
of the 2601 nodes. This reference displacement field is shown in figure ©• This target 
displacement field was differentiated by finite differences to calculate the target strains 
shown in figure a- Since the analytical expressions for the displacements are available, 
the target strains can be differentiated exactly rather than using finite differences. There 
is no analytical expression for the noise corrupted displacement fields, however, and 
therefore we have chosen to treat both noisy and reference fields in the identical manner. 
These target fields will be used as a reference to benchmark the accuracy of the fields 
reconstructed with SPREME. 
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Figure 2: Uniaxial tension test: The exact analytical solution for far-field uniaxial ten¬ 
sion applied to an infinite sheet with a bonded circular inclusion is used for verification of 
the computational implementation. In the figure, ROI indicates the Region Of Interest 
over which displacement fields are assumed to have been “measured”. 



x 

Lateral (horizontal) displacement 



Axial (vertical) displacement 


Figure 3: Uniaxial tension test. Upper row: Exact noiseless displacement field for 
uniaxial tension experiment. These represent the target displacement fields. Lower row: 
Noise-corrupted displacement fields used for verification. The noise level in the lateral 
displacements is similar to that seen in physically measured displacement fields shown 
in figure El 


Our problem assumptions include that we are given accurate axial displacements 
and noisy lateral displacements. Therefore, Gaussian noise was added to the lateral 
displacement (u x ) in order to simulate ultrasound displacement measurements, and the 
axial displacement (u y ) was left unchanged. The corrupted displacement field is shown 
in figure ([5]). 

The corrupted displacement field was differentiated by finite differences to obtain the 
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Figure 4: Uniaxial tension test results. Top row: Target noiseless strain fields from 
uniaxial tension test. Middle row: Strain fields calculated from noise-corrupted input 
displacement fields. Bottom row: Strain fields recovered from SPREME. 


strains displayed in figure Q. The e yy strain component was unchanged. The features 
in the other strain components, however, were completely obscured by noise. 

The corrupted displacement field was processed with SPREME using the following 
parameter values: a Q = lO -5 , /? = 1, 6 = 10 -8 , n = 0.5, T xx = 10 -9 , T yy — 10 4 . 
These parameters were chosen by trial and error from extensive computational testing. 
This same trial and error method was used to choose the best parameters for all the 
other computational experiments performed. Six iterations were performed to allow the 
strains to converge. The strains from SPREME are displayed in the bottom row of figure 
0- The reconstructed strains are very smooth because of the regularization on the 
momentum equation. This regularization term is very similar to a Total Variation (TV) 
regularization which tends to penalize reconstructions with high oscillations, but allows 
reconstructions to have jumps that need to be present to fit the data [Vogel, 2002| . 

Accuracy of reconstructed strains and displacements may be seen qualitatively in line 
plots passing through the center of the inclusion region, figure 0 - These plots show 
excellent agreement between the target strains and the reconstructed strains (referred 
to as SPREME Strain in the legend) in every component of strain. We recall that the 
SPREME formulation has independent variables for both displacements and strains, and 
therefore we also use this figure to evaluate consistency of the SPREME displacements 
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and strains. To see this, we use a finite difference approximation of the derivative of the 
reconstructed displacements to obtain another estimate of computed strains; these are 
labeled SPREME Disp in the legend, and are also in excellent agreement with the other 
quantities. 





Figure 5: Uniaxial tension: Line plots through the axial, lateral, and shear strain fields. 
The SPREME method produces independent estimates of displacement and strain. Thus 
each graph shows three curves. One is the target strain. The second is the strain 
output from SPREME labeled “SPREME Strain.” The third is the derivative of the output 
SPREME displacement, labeled “SPREME Disp”. Comparing the latter two provides a 
test for consistency between these fields. 


Overall quantitative accuracy of the displacements (strains) computed with SPREME 
were evaluated for both individual components and for the total vector (tensor). The 
component error is defined as: 


Strain error 



(49) 


TV \\ U \~ U i\\ / rr\\ 

Disp error = —p— tt .— (50) 

IKII 

Here u\, and e\- are the target displacement and strain components, respectively, and 
the i's and j's go from 1-2 with 1 representing the x direction, and 2 representing the y 
direction. The counter a runs over all nodes in the mesh. 

The norms calculated with these formulas were plotted as a function of iteration 
number in figure Q. For iteration zero, u was the displacement field corrupted by 
noise, and e was initialized to zero strain. The plots in the figure demonstrate that the 
errors reduce significantly with iteration number. The exact values of the errors in the 
field variables after the sixth iteration are displayed in table 0. From the table, we see 
that the error in u x (resp/. e xx ) reduced from about 50% (resp. 500%) to about 2%. 

The total error in strain and displacement were computed with the following formu¬ 
las: 


Total strain error 


Total disp error 


/ I f err 1 
/ 1 l^xx 1 

l 2 + l 

\ f err 

\ yy 

|| 2 H 

\ Ik* 1 

y W^xxi 

l 2 + 

U I 
11 yy 1 

| 2 4 

\K rr \ 

l 2 + 

1 \n.err 112 

\\ a y II 

V \\<\ 

l 2 + 

IKII 

12 


f err 112 
c xy I1 


ct 112 

xy II 


(51) 

(52) 


Here u^ rr = U{ — u\ and e\ r - r = — e\-. A plot of these total errors as a function of 

iteration number is shown in figure Q, which shows that the total errors decrease with 
iteration. 
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Displacement error plot 


Strain error plot 




Figure 6: Uniaxial tension test: Error in each component of the field variables as a 
function of iteration number. 


Strain and displacement error 



Figure 7: Uniaxial tension test: Total error in the field variables as a function of iteration 
number. 


Field variable 

Error in input field (%) 

Error in reconstructed field (%) 

u x 

48.8 

2.09 

u y 

0 

0.072 

M 

22.6 

0.759 

txx 

511 

1.78 

e yy 

0 

0.526 

t X y 

1883 

13.5 

6 

277 

1.49 


Table 1: Uniaxial tension test: Errors in the recovered field variables after the sixth 
iteration, and initial errors in the input field variables. 

3.3 Biaxial tension test 

For the second verification test, we consider the same body, this time under biaxial 
tension; c.f. figure The stress at infinity was a xx = a yy = a°° = 6 x 1CT 2 . Next, 
Gaussian noise was added to the lateral displacement field, and reconstructions were 
performed using the following parameters: a — 10 -5 , /3 = 10, <5 = 10 -8 , n — 0.5, T xx — 
50, T yy = 5 x 10 4 . These are approximately the same parameters as the last experiment, 
except with T xx significantly increased to accommodate the applied horizontal stress. 
The target and reconstructed strains at the sixth iteration are shown in figure ©• 
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Figure 8: Biaxial tension test: A second verification experiment uses the exact analytical 
solution for far-held biaxial tension applied to an infinite sheet with a bonded circular 
inclusion. 



u U.O -u.o u u.o -u.o U 

XXX 


e xx Lateral strain e yy Axial strain e xy shear strain 

Figure 9: Biaxial tension test results. Upper row: Target noiseless strain fields from 
uniaxial tension test. Lower row: Strain fields recovered from SPREME. Relatively large 
weights on noisy lateral displacement increase the error in the recovered strains relative 
to those in Figure [4j 


The features in the lateral and shear strains are recovered, but we see significantly more 
error in this test case than in figure s- That error is quantified in table <©• The table 
shows clearly that the errors in the lateral displacement, and lateral and shear strains 
reduced significantly from the errors in the input held. This means that the method is 
able to work also for this type of deformation. The method did not work as well as it 
did for the case when the domain was under uniaxial tension, however. This is most 
likely due to the need to use relatively large weights on the noisy lateral displacements 
to capture the effect of signihcant lateral normal stress applied at inhnity. Thus was the 
effect of the noise in that displacement component amplihed. Nevertheless, the error in 
e xx reduced from about 500% to about 15%, a very signihcant reduction. 
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Field variable 

Error in input field (%) 

Error in reconstructed field (%) 

u x 

51.2 

6.63 

u y 

0 

0.103 

M 

36.6 

3.83 

^xx 

495 

15.4 

e yy 

0 

1.04 

e xy 

1370 

59.3 

e 

422 

13.9 


Table 2: Biaxial tension test: Errors in the recovered field variables after the sixth 
iteration, and in the input field variables for the biaxial tension experiment with noise. 


4 Validation 


In the verification section, we demonstrated that SPREME correctly reconstructs plane 
stress displacement fields. Here in the validation section, we demonstrate that the 
SPREME implementation of the plane stress model applies to physical systems of interest. 
The displacement data used for the test was measured within a gelatin phantom under 
uniaxial compression as described in |Dor d et al., 2012 . This displacement data was an¬ 


alyzed previously by an iterative elastic modulus inversion algorithm |Dord et al., 2012| . 
The results from that study, and a smoothed version of the measured displacements, will 
be used as a reference to compare to the results produced by SPREME. 


4.1 Phantom tests 

The displacement data used for the test was measured within a gelatin phantom un¬ 
der uniaxial compression as described in |Dord et al ., 2012 . Here, for completeness of 
presentation, we briefly describe the data collection procedure. An agar and gelatin 
phantom was manufactured of to have the acoustic and mechanical properties of soft 
tissue using methods described in |Pavan et a l., 2010 . The phantom was a 100mm cube 
containing four spherical inclusions of diameter 10mm. The inclusion centers are copla- 
nar in a horizontal plane separated by 30mm center to center distance. Each of the 
inclusion has a different material Young’s modulus [Dord et al., 2012| . 

4.1.1 Ultrasound imaging and displacement estimation 

The phantom was imaged with a Siemens SONOLINE Antares (Siemens Medical So¬ 
lutions USA, Inc, Malvern, PA) clinical ultrasound system, with a linear ultrasound 
transducer array (Siemens VFX9-4) which was pulsed at 8.89MHz [Dor d et al., 20121. 
A 15cm x 15cm compression plate, much larger than the phantom surface, was attached 
to the ultrasound transducer to help generate approximately uniaxial deformation of the 
phantom. The phantom was first imaged before any deformation was applied, and it 
was imaged after every 0.5% strain up until a total strain of 20% with respect to the 
phantom’s height |Dord et al., 2012] . During the imaging, Radio Frequency (RF) data, 
representing the spatial distribution of the backscattered pressure field, was recorded. 
A modified block matching algorithm was used to estimate the displacement field from 
the RF data | Jiang and Hall, 201 1| . Because the models used in this paper depend on 
the small strain assumption, we use only the measured deformation fields corresponding 
to small strain, specifically, we choose 1.5% overall strain. 
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Figure 10: Axial displacement fields. First row: measured displacement field, second 
row: displacements reconstructed by SPREME, third row: reference displacements from 
an iterative inversion algorithm. 


4.1.2 Reference results 

We compare our reconstructed displacement fields to two reference fields. 

The first reference displacement is simply the spatially smoothed version of the mea¬ 
sured lateral displacement field. This smoothing is achieved by performing a local spatial 
average with corrections at the boundaries. A smoothed displacement field reveals its 
main features, but, of course, obscures detail. 

A second reference displacement field is computed as the result of a second inverse 
problem to reconstruct the linear elastic shear modulus, as described in |Goenezen et ah, 2012 
To that end, given a current guess of the modulus field, a forward elasticity problem 
is solved driven by assumed boundary conditions. The modulus field is updated to 
minimize the difference between measured and predicted displacements. Only axial dis¬ 
placements are used in the minimization. 

4.1.3 Phantom results and discussion 

The displacement for each target was downsampled to a 63 x 54 grid, which is the same 
grid that was used to generate the reference results |Dord et ah, 2012| . Images of the 
measured, smoothed, reconstructed, and reference displacements are shown in figures 
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Figure 11: Lateral displacement fields. First row: measured displacements; second row: 
a smoothed version of the measured displacements; third row: displacements recon¬ 
structed by SPREME; fourth row: reference results from an iterative inversion algorithm. 
We note general agreement in gross features between all rows for targets 1 and 2. For 
targets 3 and 4, there is disagreement between the reference and the other fields, which 
we attribute to inaccuracies in assumed boundary conditions used to create the reference 
solution. 


(10) and 


md {lT| . Th 

— irFT o _ 


The parameters used for SPREME processing were: T xx — 10 
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± yy 


1, a — 10 /? = 10, S = 10 , n = 0.5. Eleven iterations were performed for each 

th 

target, though it was observed that the strains typically converged by the 5 iteration. 

We note that the reconstructed axial displacements match the measured and refer¬ 
ence axial displacements very closely. The reconstructed lateral displacements, however, 
were less similar to the reference displacements. This difference is most noticeable in 
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targets 3 and 4. It is interesting to note that the general features in the SPREME dis¬ 
placement images for both these targets are more similar to the smoothed measurements 
than are the general features in the reference displacement images. This indicates that 
the reference lateral displacements are in error, and this is likely due to errors in the 
assumed boundary conditions used to compute the reference results. 
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Figure 12: Axial strain fields. First row: measured strains; second row: strains recon¬ 
structed by SPREME; third row: reference results from an iterative inversion algorithm. 
We note qualitative and quantitative agreement between the rows; we note also the 
presence of boundary artifacts in the reference fields that are not present in the SPREME 
fields. 


Axial strain images (e yy ) for each target are shown in figure ( [l2| ), while the lat¬ 
eral strain fields ( e xx ) are in figure (13). The measured (resp. reference) strains were 
computed by differentiating the measured (resp. reference) displacement fields using a 
finite difference approximation. Three observations may be made. First, except for the 
lateral measured strain, all strain fields show clearly the presence of a stiff inclusion of 
the same size and shape, and with similar strain contrasts. The high noise level in the 


2 We recall that the iterative inversion code performs a forward elasticity solve at every iteration. 
This solve requires assumed displacement and traction boundary conditions. The lateral direction is 
assumed to be traction free. The measured axial displacements were fixed as boundary conditions on 
all four edges of the phantom. 
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Figure 13: Lateral strain fields. First row: measured strains; second row: strains recon¬ 
structed by SPREME; third row: reference results from an iterative inversion algorithm. 
We note general agreement of the main features between the SPREME and reference re¬ 
sults; boundary noise present in reference results are not apparent in SPREME results; 
SPREME fields show higher dynamic range than reference fields. 


measured lateral displacement fields, however, is magnified through differentiation to 
a point where it obscures all features that might be present in the data. Second, we 
note boundary artifacts in the reference strain fields; these are caused by the assumed 
boundary conditions in the reference reconstruction method. The SPREME strain fields 
are relatively free of such artifacts. Third, strain contrast or dynamic range is slightly 
(resp. significantly) higher in the SPREME axial (resp. lateral) strain reconstructions than 
in the reference fields. This is perhaps due to deleterious effects of regularization used in 
the reference method, which tends to diminish modulus, and thus reconstructed strain, 
contrast. Similar observations may be made when comparing shear strains, figure (14). 
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Figure 14: Shear strain fields. First row: measured strains; second row: strains recon¬ 
structed by SPREME; third row: reference results from an iterative inversion algorithm. 
We note general agreement of the main features between the SPREME and reference re¬ 
sults; boundary noise present in reference results are not apparent in SPREME results; 
SPREME fields show higher dynamic range than reference fields. 


5 Application to in-vivo image data 

In this section, we test SPREME with in-vivo displacement data measured from patients 
with breast masses. The main goal of this test is to evaluate whether SPREME is suf¬ 
ficiently robust to provide useful results with clinical data, which is the motivating 
application. The in-vivo displacement data used for the tests has been processed with 
an iterative inversion code in a different study described in |Goenezen et ah, 20121. The 
results from that study will be used as a reference to benchmark the results produced by 
SPREME. These results contain displacement fields from 5 biopsy proven fibroadenomas 
and 5 biopsy proven invasive ductal carcinomas, which represent the most common forms 
of benign and malignant breast tumors. For the present purposes of demonstration, in 
this section we present results from one of each; the reconstructed displacement and 
strains for the other eight cases are shown in appendix ([B]). For all the cases treated, we 
work on the same mesh as used in |G oenezen et ah , 20121, and focus on displacement 
fields corresponding approximately to roughly 1% overall strain in order to stay in the 
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linear range. The input measured displacement fields used are shown in Figure (15) 


5.1 In-vivo results and discussion 
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Figure 15: Displacement fields measured in-vivo in the breast as described 
|Goenezen et ah, 2012 . Upper row: Fibroadenoma, a benign breast tumor (FA 3) Lower 
row: Invasive Ductal Carcinoma, a malignant breast tumor (IDC 2). 


The parameters used for SPREME processing of the clinical data were: T xx — 10 -3 , T yy 
1, ck = 5x 10 -4 , (3 = 10, S = 10 -8 , n = 0.5. The reconstructed displacement and strain 
fields obtained after the eleventh iteration are shown in figures (16) and along with 
the reference fields. As with the validation study, we compare to two reference fields: 
the smoothed measured lateral displacement fields, and the results obtained by an iter¬ 
ative inversion algorithm published in |Goenezen et ah, 2012 . Generally speaking, we 
see excellent agreement between all three fields in the axial displacement and strain 
components. 

In the lateral displacement and strain fields, however, there is less consistency. The 
SPREME lateral displacement fields tend to reproduce the gross features seen in the 
smoothed lateral displacement measurements; the reference results, however, sometimes 
agree (c.f. figure ©), and sometimes disagree (c.f. figure (©) with the other two. 
The reference displacements seem to indicate an outward motion at the left and right 
side of the imaged domain in all the tumor cases (even in the 8 other cases shown 
in the appendix). The smoothed measurements and the predicted displacements from 
SPREME, however, do not indicate this type of motion in all the tumor cases. This 
outward motion predicted by an iterative inversion code is probably an artifact due to 
the assumed boundary condition of traction free sides on the left and right sides of the 
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Figure 16: Displacement and strain images for fibroadenoma 3. We note that the 
reconstructed SPREME lateral displacements reproduce the gross features seen in the 
smoothed measurements, while the reference method does not. 


imaged domain, and an implied symmetry. These assumptions may not be always valid 
in practice. 

In contrast to the lateral displacement fields, the lateral strain fields recovered from 
SPREME were remarkably similar to the reference strains in all cases. This is surprising 
because the reconstructed lateral displacements were not always similar to the reference 
lateral displacements. This implies that the dominant difference between the SPREME 
and reference lateral displacement may be a rigid body mode. 

6 Summary and conclusions 

In this paper we consider the problem of reconstructing both components of a 2D vector 
displacement field in a heterogenous elastic medium, given a precise measurement of 
one of the two components and a noisy measurement of the other. This problem is 
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Figure 17: Displacement and strain images for invasive ductal carcinoma 2. In this ex¬ 
ample, we see more consensus between the gross features of the smoothed measurements, 
the reconstructed SPREME lateral displacements, and the reference lateral displacements. 


motivated by applications in quasistatic ultrasound elastography, in which technological 
limitations allow the very precise measurement of one component of displacement, but 
impede the measurement of others. 

Under the assumption that the material is piecewise homogenous, we proved that the 
momentum equations and knowledge of one component of the displacement field deter¬ 
mines the other displacement component up to four undetermined coefficients. We then 
went on to derive a variational formulation that can exploit the condition of piecewise 
homogeneity without other a priori knowledge of modulus distribution. 

We then showed verification, validation, and application of a finite element implemen¬ 
tation from this variational formulation. The verification was performed on simulated 
data and showed that the method described here converges quickly, and dramatically re¬ 
duces error in the noisy (uncertain) displacement component. This came at the expense 
of a slight increase in noise in the precise components. Nevertheless, in the verification 
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studies conducted, the overall strain error reduced from roughly 0(1) (several hundred 
percent) to 0( 10 -1 — 10 -2 ) (a few percent). The validation studies showed that this 
method produces realistic displacement fields from measurements in tissue mimicking 
phantoms. The predicted displacements tend to be more consistent with measurements 
than those of a reference method in the literature. Finally, the application of SPREME 
to data collected in-vivo demonstrates that SPREME is sufficiently robust to work in the 
application domain that originally motivated its development. 

Further improvements in lateral displacement estimation with SPREME are possible. 
Since SPREME is a post-processing method, it can be used in conjunction with the 
other lateral displacement estimation techniques outlined in the introduction to obtain 
more precise estimates of the displacement field. Since these methods produce improved 
measurements of the lateral displacements, then relatively large weights can potentially 
be used on the lateral displacements when processing them with SPREME. This processing 
can be done efficiently because SPREME converges in very few iterations. 

An area that needs further work is finding a systematic method to identify good 
choices of the algorithmic parameters in SPREME. Regularization parameter selection 
is an open and active area of research in inverse problems, and is not addressed in the 
present study. 

The SPREME approach may be generalized to other applications, as described in 
|Babaniyi, 2012 . The idea is that one identifies the form of the momentum equation that 
results for a homogeneous material property distribution. This equation is then enforced 
almost everywhere, as indicated in equation ©• See [Babaniyi, 2012| for details. 
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A Formulation without compatibility 


(k — l ) , 

Given u m (x,y), x £ Q and e , find e k that minimizes: 


ir[e \=tt 0 +tt r 


(53) 


where: 


*° = 2 
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and: 
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(V-A (e fc )) ; 


(V- A(e (fc - 1) )) 2 + <5 


dfl. 


(55) 


In equation (54), T is a fourth order weighting tensor that allows more importance 


to be placed on the more accurate strain component. So in ultrasound measurements 
where the axial strains (e yy ) are typically more accurate than the lateral and shear 
strains (e xx ,e xy ), then the T yyyy weights will be larger than T xxxx and T a 


xyxy 


Equation ([53]) was minimized with respect to e to get the weak form using the same 
method outlined in section 12.21 The weak form was then discretized with finite element 
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Figure 18: Reconstructed strains at fifth iteration. 


bilinear shape functions to get a linear system of equations. A finite element within an 
in-house FEM code was created to solve the linear system of equations. 

This formulation was tested with the perfect input displacements shown in figure 

The following parameters were used to run the code: Tyyyy - 1, T XXXX - T X y X y - 

1 x 10 -9 , a = 1 x 10 -6 , 5 = lx 10 -8 , n = 0.5. These weights were chosen in order 
to simulate ultrasound measurements where the lateral and shear strains are not well 
known. Five iterations were performed, and the reconstructions obtained at the fifth 
iteration are shown in figure (18). The lateral and shear strains look very different from 
the target distributions shown in figure 0- 


Field variable 

Error in reconstructed field (%) 

^xx 

125 

e yy 

1.41 

e xy 

695 


Table 3: Errors in the recovered strains after the fifth iteration. 


To quantitatively determine how accurate the reconstructions were, the errors be¬ 
tween the reconstructed strains at the 5th iteration and target strains were computed 


using the formula shown in equation (49). The results of this computation are displayed 
in table Q. These results show that the lateral and shear strains are not accurate. 

Another method used to check the output strains from the FE code is the compati¬ 
bility equations which are shown below Atkin and Fox, 2005 : 


£xx,yy T £yy,xx — 2 £xy,xy 


(56) 


The equation used to check the compatibility condition is slightly different from the 


equation shown in (56). The derivation of the equation used is shown below: 


2 £xy — ^x,y T ^y,x 

(57) 

2 £xy ^x,y — 

(58) 

{^^xy ^x,y\y ~ ^ y,xy 

(59) 

^ x,y\y ^ yy,x = V’ 

(60) 

a measure of incompatibility. 

Setting rj — 0 in (60), 


and taking its derivative with respect to x, we see that it is the same as equation (|56|). 
Therefore by knowing u 


e xy , and e yy , we can use equation (60) to check that the strains 


from the FE code satisfy the compatibility equation, r] was evaluated pointwise within 
the domain for each of the experiments performed. To get a sense of how well the 
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Table 4 : Compatibility results using various weights. 


compatibility equation is satisfied for each experiment, we calculate the L 2 norm. The 
formula used to evaluate the L 2 norm is defined below: 



The L 2 norm of 77 was computed for a large range of experiments where the weights 
were varied, and the parameters ct, 5 , and n were fixed to 1 x 10 -6 , 1 x 10 -8 , and 
0 . 5 , respectively. The results from these experiments are shown in table 0- The re¬ 
sults on the table demonstrate that the compatibility equation becomes increasingly 
violated as T xxxx and T xyxy are reduced, and the level of violation is maximum when 
T XX xx — T xyX y = 10 -5 . The mesh size for the experiment is 0.02 therefore any value of 
| lr /11 smaller than this mesh size can be neglected as discretization error. When T xxxx 
and T xyxy becomes smaller than 10 -3 , then the compatibility equation becomes vio¬ 
lated for the experiment. This means that this alternate formulation cannot be used to 
calculate correct strains when the weights for the lateral and shear strain components 
are too small. When the weights for some of the strain components are small, the FE 
code constrains the corresponding strains with small weights to satisfy the equilibrium 
equations, but the strains are not forced to come from a single valued continuous dis¬ 
placement field. This formulation will therefore not work on real measured data where 
we would wish to use small weights to exclude the strain components calculated from 
the imprecisely measured displacement component. 

By way of comparison, the SPREME formulation yields ||r/|| ~ 9 x 10 -3 . 


B Clinical data reconstructions 

In this section, for completeness, we show the reconstructions obtained from processing 


the rest of the in-vivo data from G oenez en et ah, 2 012| . The observations made in 
earlier apply to these results as well. 
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Figure 19: 


Displacement and strain images for fibroadenoma 1 
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Figure 20: Displacement and strain images for fibroadenoma 2 
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Figure 21: Displacement and strain images for fibroadenoma 4 
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Figure 22: Displacement and strain images for fibroadenoma 5 


30 























































Lateral strain Axial strain Axial disp. Lateral disp. 

y (mm) y (mm) y (mm) y (mm) 


Smoothed measured lateral disp. (u ) 



0 10 20 30 

x (mm) 



0 10 20 30 

x (mm) 


Measured axial strain (e ) 

V yy > 



x (mm) 


Measured lateral strain (e ) 

xx x 10" 



0 10 20 30 

x (mm) 


Measured 


Reconstructed lateral disp. (u ) 


Iterative lateral disp. (u ) 



Reconstructed axial disp. (u ) 


Iterative axial disp. (u ) 



Reconstructed axial strain (e ) 



Iterative axial strain (e ) 



Iterative lateral strain (e ) 



Reconstructed (spreme) 


Reference 


Figure 23: Displacement and strain images for invasive ductal carcinoma 1 
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Figure 24: Displacement and strain images for invasive ductal carcinoma 3 
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Figure 25: Displacement and strain images for invasive ductal carcinoma 4 
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Figure 26: Displacement and strain images for invasive ductal carcinoma 5 
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